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We present an approach for generating local numerical basis sets of improving accuracy for first-principles 
nanoplasmonics simulations within time-dependent density functional theory. The method is demonstrated 
for copper, silver, and gold nanoparticles that are of experimental interest but computationally demand¬ 
ing due to the semi-core d-electrons that affect their plasmonic response. The basis sets are constructed 
by augmenting numerical atomic orbital basis sets by truncated Gaussian-type orbitals generated by the 
completeness-optimization scheme, which is applied to the photoabsorption spectra of homoatomic metal 
atom dimers. We obtain basis sets of improving accuracy up to the complete basis set limit and demonstrate 
that the performance of the basis sets transfers to simulations of larger nanoparticles and nanoalloys as well 
as to calculations with various exchange-correlation functionals. This work promotes the use of the local basis 
set approach of controllable accuracy in first-principles nanoplasmonics simulations and beyond. 


I. INTRODUCTION 

Plasmonics attracts increasing interest due to its tech¬ 
nological relevance in numerous applications, such as bio¬ 
chemical sensing, 1 sub-wavelength light manipulation, 2 
and photovoltaics. 3 Plasmon resonances in metal nano¬ 
particles can be qualitatively understood by classi¬ 
cal electromagnetism, but the accurate description of 
nanometer-size particles or systems with features in the 
subnanometer range requires more elaborate approaches. 
In these systems the plasmonic response is affected by 
quantum effects, such as electron spill-out at the surface 
and electron tunneling. 4 A number of recent studies has 
demonstrated that the regime where these effects start 
to play a role is experimentally accessible. 5-9 

To study quantum effects in the plasmonic re¬ 
sponse computationally, one often resorts to time- 
dependent density functional theory 10 (TDDFT) simu¬ 
lations. Qualitative understanding on the quantum ef¬ 
fects in nanostructures can be obtained within the jel- 
lium approximation, 11-14 but it cannot capture the im¬ 
portant atomic structure effects. 15 In addition, the jel- 
lium model only describes simple metals such as sodium, 
where the optical response is determined by the valence 
s-electrons. However, the experimentally relevant mate¬ 
rials for plasmonics are usually coinage metals; copper 
(Cu), silver (Ag), and gold (Au). 5-9,16,17 In these metals, 
in addition to the outermost s-electrons, also the semi¬ 
core d-electrons participate in the response. Although 
the effects due to the d-electrons can be accounted for 
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in an approximate manner in the jellium model, 18 first- 
principles models are necessary to obtain accurate re¬ 
sults. 

Coinage metal systems have been studied through nu¬ 
merous TDDFT simulations, including metal nanopar¬ 
ticles of different shapes, 19-26 nanoalloys, 27-31 protected 
metal clusters, 32-36 and nanoparticle dimers. 37 However, 
TDDFT simulations for these systems are computation¬ 
ally demanding. Even though the calculations can be 
speeded up with the frozen-core approximation, coinage 
metals require explicit calculation of the semi-core d- 
electrons in addition to the s-electrons, resulting in 11 
electrons per atom in calculations, in contrast to, e.g., 
sodium where it usually suffices to treat only the single 
3s-electron per atom. Consequently, simulated systems 
have typically been restricted to the maximum size range 
of 100-200 coinage metal atoms, 19-23,27-33,35-37 with a 
few studies presenting larger systems, such as a Au 263 
nanorod, 26 a Ag 272 nanoshell, 24 and a thiolate-protected 
Au 3 i 4 . 34 The calculations have employed either real- 
space grid codes 26-28,33-36 or the linear combination of 
atomic orbitals (LCAO) approach. 19-24 ’ 29-32,3 ' Recently, 
a new LCAO-TDDFT implementation was developed, 25 
allowing to push the accessible system size close to the 
classical limit (Ag 561 presented in Ref. 25). 

A serious problem of the LCAO approach is that it 
is prone to errors due to basis set incompleteness — a 
problem which can be straightforwardly tackled in the 
real-space grid methods. Nevertheless, this issue has not 
been extensively discussed in previous nanoplasmonics 
studies using the LCAO approach. Instead, a reasonable 
accuracy of the results has been checked by calculations 
of test systems with larger basis sets of the available ba¬ 
sis set series 19,20,24 or by comparing to real-space grid 
results. 25 
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The basis set issue is complicated by the fact that con¬ 
ventional basis sets are typically optimized for ground- 
state energy calculations. 38 For other properties, such as 
dipole moments, excited states, or plasmonic charge den¬ 
sity oscillations, these basis set series are not expected 
to yield quickly converging results. 39-43 In the case of 
photoabsorption spectra, the accurate description of the 
dipole moment of the excited states is essential. This 
necessitates inclusion of diffuse functions (i.e., functions 
with large spatial extent) in the basis set. Diffuse func¬ 
tions are not present in the energy-optimized basis sets 
because of their minor contribution to the ground-state 
energy of electrically neutral systems, but they can be 
generated into these basis sets by, e.g., minimizing the 
energy of anions. 38 Alternatively, unoccupied atomic or¬ 
bitals can be included in the basis set, which has been 
found to improve the description of the photoabsorption 
of metal nanoparticles. 25 ’ 44 However, these approaches 
are not guaranteed to be optimal for extending the ba¬ 
sis sets beyond conventional ground-state energy calcu¬ 
lations. 

In the present work, we show that efficient basis sets 
specifically optimized for describing the plasmonic opti¬ 
cal response can be systematically generated using the 
completeness-optimization (CO) approach. 41 CO is a 
black-box procedure for generating basis sets for any 
property at any level of theory. It has been previ¬ 
ously used to generate all-electron Gaussian-type or¬ 
bital (GTO) basis sets for calculating magnetic 41 ’ 45-50 
and magneto-optic° 1-57 properties as well as the elec¬ 
tron momentum density. 42,58 In this work, we present a 
straightforward extension of the CO formalism to semi- 
numerical basis sets by combining numerical atomic or¬ 
bitals (NAOs) with truncated numerical Gaussian-type 
orbitals (NGTOs). The NGTOs are selected by a re¬ 
cently developed automatic CO procedure 42,43 to aug¬ 
ment NAO basis sets with the necessary diffuse and po¬ 
larization functions. We demonstrate the applicability 
of the scheme for describing collective plasmonic exci¬ 
tations in coinage metal clusters. We optimize the ba¬ 
sis sets to reproduce the photoabsorption spectra of ho- 
moatomic dimers, and show that the generated basis sets 
are transferable to larger nanoparticles and to different 
chemical environments in nanoalloys, as well as to differ¬ 
ent exchange-correlation functionals. 

The paper is organized as follows. In Sec. II, we give an 
overview of the used methodologies — TDDFT, LCAO, 
and CO. In Sec. Ill, we describe our implementation, 
and in Sec. IV we demonstrate the performance of the 
basis set generation and test the transferability of the 
generated basis sets. We conclude the study in Sec. V. 

II. METHODS 

A. Time-dependent density functional theory 

TDDFT is a well-established formulation of the time- 
dependent many-body Schrodinger equation in terms of 


the time-dependent electron density. 10 The theory is usu¬ 
ally applied within the Kohn-Sham (KS) description of 
density functional theory (DFT), 59,60 which models the 
interacting many-electron system as a non-interacting 
system in an effective potential. In this approach, the 
complicated many-body interactions are described by the 
so-called exchange-correlation (xc) functional. The time- 
dependent xc functionals are usually treated in the adia¬ 
batic limit, i.e., an instantaneous time-dependent density 
is used as input for the ground-state functional. 61 

The dynamical response, and in particular, the pho¬ 
toabsorption spectrum of a given system, can be cal¬ 
culated in two formally different, but equivalent man¬ 
ners within the TDDFT framework. First, it can be 
obtained from the time-dependent dipole moment that 
is recorded during the explicit real-time propagation of 
the KS-orbitals that have been excited from the ground 
state by a (5-pulse perturbation. 62 Second, the excitations 
of the system can be calculated by formulating the lin¬ 
ear density response to an external perturbation in the 
frequency space, yielding the Casida matrix equation. 63 

In this work, we use both the time-propagation and 
Casida schemes for calculating photoabsorption spec¬ 
tra. We employ the open source GPAW program 64-69 
in TDDFT calculations. GPAW uses the projector aug¬ 
mented wave (PAW) method' 0 for freezing the inert core 
electrons and for obtaining smooth pseudo-wave func¬ 
tions in the vicinity of the nuclei. The simulations ex¬ 
plicitly include only the outermost electrons, he., for the 
coinage metals, the semi-core d-electrons and the va¬ 
lence s-electrons (11 electrons per atom in total). The 
element-specific PAW transformations are constructed at 
the scalar-relativistic level of theory. Thus, relativis¬ 
tic effects, especially important for gold,' 1 are included 
implicitly in the calculations through the PAW trans¬ 
formation. For the present study, GPAW has the ad¬ 
vantage that it can describe wave functions either on a 
real-space grid 64,65 or within the LCAO approach. 25,66 
In both modes, uniform real-space grids are used for rep¬ 
resenting electron densities and potentials. These two 
modes of operation share a significant portion of compu¬ 
tational framework within the program, which allows us 
to compute grid-based spectra and LCAO spectra with 
minimal sources of differences apart from the representa¬ 
tion of the wave functions. 


B. Linear combinations of atomic orbitals 

In the LCAO approach the single-electron KS wave 
function is expressed as 

atoms N a 

w=EE c “ itf>, (!) 

a v—l 

where c“ are the expansion coefficients for basis functions 
\\v) centered on the atom a, N a denoting the amount of 
basis functions on that atom. In a coordinate system 
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centered on the atom a, an associated basis function is 
written as a product of a radial function r n i(r) and a 
spherical harmonic Yi m (9,ip), 

{r\xl) = X*(r) = ^iA r ) Y ^rnA e ^)- ( 2 ) 

Above, v is a symbolic index over the combinations of 
n„, and m u . In this work, radial functions are taken 
to be either NAOs or NGTOs as described in Sec. III. 

The main advantage of the LCAO approach is that a 
sufficiently accurate description of the wave function can 
often be achieved with a small number of basis functions. 
In addition, in the case of truncated basis functions, the 
number of overlapping basis functions is usually small, 
which enables efficient computations. The main draw¬ 
back of the LCAO approach is that it is prone to errors 
due to the incompleteness of the basis set. Thus, it is 
important to use basis sets that are flexible enough for 
describing wave functions accurately in the regions that 
are essential for the property in question. 

C. Completeness-optimization 

CO is a general approach for generating optimal basis 
sets for any chosen property . 41 The method is based on 
the concept of the completeness profile 40 that is defined 
as 


to find the optimal limits of the intervals, while keeping 
Yi(a ) « 1 within the intervals. 

In this work, we use Gaussian basis functions char¬ 
acterized by their exponents a and employ the auto¬ 
matic CO procedure 42 ’ 43 implemented in the open source 
ERKALE program .' 2,73 The procedure is based on sys¬ 
tematic trial and error searches for determining the com¬ 
plete basis set (CBS) limit of the property, as well as the 
CBS itself. The algorithm is divided in two phases. First, 
in the extension phase, the CBS is found by progressively 
extending the basis set with the most important functions 
until the property is converged. Each addition of a ba¬ 
sis function corresponds to an extension of the intervals 
[a[ nln , of ax ] or decrease of the tolerance 77 . Second, in 
the reduction phase, the least important basis functions 
are repeatedly pruned from the basis, shrinking the 07 
intervals or increasing 77 . During the optimization, the 
relative importance of the addition or removal of a basis 
function is determined by a user-defined error metric. In 
every step of the algorithm, the exponents of the Gaus¬ 
sian basis functions are determined by maximizing the 
completeness profile Y) within the current ct; interval. As 
a result of the reduction phase, a systematic sequence of 
basis sets of decreasing size and accuracy is obtained. 

III. IMPLEMENTATION 


N 

Y i{ot)= E {9i(<x)\Xn) S~^ {x v \gi{a)), (3) 
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where \xv) are the basis functions, S ~ 4 is the {p,v) 
element of the inverse overlap matrix orthonormal- 
izing the basis, and gi(a) is a primitive test func¬ 
tion usually taken to be a normalized Gaussian prim¬ 
itive gi{a) oc r l e~ ar Yi m (d,ip). For basis functions of 
the form of Eq. (2), the inner product {gi(ct)\xv) oc 
(: r lo, e~ ar \<j>n„ij) Thus, each value of the an¬ 

gular momentum l yields a different completeness profile 
Yi(a), whereas all m values of a given l yield the same 
profile. 

The profile essentially measures the validity of the res¬ 
olution of the identity operator 

N 

E ( 4 ) 

fi^—l 

which would be exact for a complete basis set. 
Correspondingly, the completeness profile satisfies 
0 < Yi(a) < 1. 

The idea in CO is to generate a basis set that has 
Yi(a) ~ 1 within the intervals [a“ m , a™ ax ] that are im¬ 
portant for the property in question. This is accom¬ 
plished by optimizing the parameters of the basis func¬ 
tions so that Y, is maximized within the intervals. The 
number of basis functions needed for this depends on the 
tolerance 77 for deviations from unity in Y] within the 
intervals . 41-43 The task of a practical CO algorithm is 


A. Application of the completeness-optimization in 
nanoplasmonics 


The CO routine 43 has a general interface for basis set 
generation. We have written a wrapper program that 
implements this interface and calls the GPAW program 68 
for photoabsorption spectrum calculations. The practical 
workflow between the programs is as follows. The opti¬ 
mization routine forms trial CTO basis sets and feeds 
them to the wrapper. The wrapper transforms the GTO 
sets into NGTOs, prepares the input for GPAW, and sub¬ 
mits the GPAW calculations. Once the GPAW jobs have 
completed, the wrapper reads in the results and returns 
them to the optimization routine, which interprets them 
through the error metric and uses the information to up¬ 
date the basis set and generate new trials. 

The error metric for determining the effect of modi¬ 
fying the basis set is defined as follows. The error e of 
a photoabsorption cross-section t 7 a b s (w) from the corre¬ 
sponding reference spectrum cr^ f s (w) (given by an earlier 
CO step) is defined as 


e = min 
s 


Tibs (- 0 ' 


S) — (7 


ref 

abs 




1/2 


+e W7 ) 2 _ l 


(5) 


where w m in and w max are cut-off energies, (e^/ 7 -* — 1 ) is 
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a penalty term for a constant energy shift S in the spec¬ 
trum, and 7 describes the stiffness of the penalty func¬ 
tion. The error metric also depends implicitly on the pa¬ 
rameters used to broaden the discrete TDDFT spectrum 
to model the finite lifetime of excitations as well as tem¬ 
perature effects and instrumental resolution. Without 
explicitly allowing an energy shift <5 in the definition, the 
error measure would be much more sensitive to offsets in 
energy than to changes in intensity. The energy shift <5 is 
penalized through the penalty term, in which the param¬ 
eter 7 provides a sliding scale to balance the sensitivity of 
the measure between energy and intensity. With 7 —> 0 + 
the penalty term becomes extremely stiff and S = 0 al¬ 
ways minimizes the error metric, in which case the met¬ 
ric becomes the usual normalized L 2 measure. The error 
metric is also used for benchmarking the quality of the 
obtained basis sets. In this case, the reference spectrum 
is obtained from a real-space grid calculation. However, 
it is emphasized that during the optimization the grid 
reference is not employed by any means. 

A practical issue concerning the extension of the CO 
routine to the frozen-core approximation is that the CO 
routine was initially designed for all-electron calculations. 
In the present work, the basis sets only represent the out¬ 
ermost s- and d-electrons of the coinage metal atoms, as 
the core electrons are described implicitly. Thus, the as¬ 
sumptions made for the composition of the basis set in 
the all-electron case 43 have been relaxed here, allowing 
for a non-monotonic amount of GTOs on consecutive an¬ 
gular momentum shells during basis set reduction. 


B. NAO+NGTO basis sets 

The benefit of NAOs is that they are not restricted to 
be of any analytical form. Thus, they should excel in 
the highly structured atomic core region, which is only 
slightly affected by the surrounding chemical environ¬ 
ment. Although the PAW transformation results in a 
smoother pseudo-wave function in this region, the NAOs 
are still pre-eminent for describing the atomic ground 
state. 

While the generation of minimal-basis NAOs is 
straightforward, the case for polarization and multiple- 
£ functions is not as clear. NAOs do not hold advantages 
for these functions, as the forms of the functions are not 
generally known and the polarization changes from one 
system to the other. NAO polarization and multiple-^ 
functions are typically generated from gas-phase atomic 
wave functions by splicing the radial function ' 4 and/or by 
including bound unoccupied orbitals, but also other sys¬ 
tematic approaches for generating NAO basis sets have 
been developed . 75-78 Here, (N)GTOs have a definite ad¬ 
vantage, as systematic sets of functions of any angular 
momentum can easily be generated. In the default ba¬ 
sis sets of GPAW, NGTOs are employed in addition to 
numerical orbitals and splicing . 66 

We employ the following strategy for basis set gener¬ 


ation. NAOs are used to represent the atomic ground 
state, for which they hold the definitive advantage. This 
minimal basis is then augmented with NGTOs that are 
generated by the CO routine, adding the desired polar¬ 
ization and diffuse functions. 

The CO routine treats analytical GTOs, whereas in 
this work, numerical basis functions are used. Thus, the 
analytical GTOs are smoothly truncated to NGTOs in 
the wrapper program by a second-order polynomial so 
that the radial part becomes 

G{r) = Ar l (e~ ar2 ^ (a - br 2 )), ( 6 ) 

where A is a normalization constant, and the constants a 
and b are chosen so that G(r c ) = G'(r c ) = 0 at the cut-off 
radius r c . The cut-off radius is determined by requiring 
that the NGTO differs point-wise from the exact GTO 
oc r l e~ ar at most by 10 -3 au -3 / 2 . This strict criterion 
distorts the Gaussian shape negligibly, but, on the other 
hand, leads to large r c values, which has a detrimental 
effect on the computational cost. 

Because of the negligible difference between GTOs and 
NGTOs, the analytical form is used in the optimization of 
the Gaussian primitives, he., in the maximization of the 
completeness profile Eq. (3) determining the exponents 
of the primitives. Additionally, the underlying NAOs are 
neglected in Eq. (3). Although NAOs (and truncated 
NGTO forms) could be straightforwardly included in the 
optimization of the completeness profile (Eq. (3)), the ef¬ 
fects due to the explicit account of the NAOs at this stage 
are expected to be small due to the different asymptotic 
forms of NAOs and GTOs. The incorporation of the un¬ 
derlying NAO basis set is done in the wrapper program 
so that its presence is invisible to the CO routine and the 
major effect of NAOs rises implicitly through the spec¬ 
trum calculation. 

The following ad hoc restrictions are imposed for the 
allowed NGTOs. At least 90% of the norm of the function 
must be within a sphere of radius 6 A around the nucleus 
and at most 90% of the norm may be within 0.6 A. 80 
The first condition results in the rejection of extremely 
diffuse functions that would require impractically large 
simulation grids, and which would cause severe numeri¬ 
cal problems due to linear dependencies in extended sys¬ 
tems. The second condition ensures that even the tight¬ 
est functions can be faithfully mapped to the real-space 
grid used to describe their contributions to the electron 
density. However, tight functions are not usually needed 
anyhow because of the PAW transformation. 


C. Numerical parameters 

The CO was started from an initial NAO+NGTO ba¬ 
sis set composed of three radial functions on the s-, p-, 
and d-shells, totaling 27 basis functions per atom. The 
initial GTO basis set was energy-optimized 43 for the gas- 
phase atom in question. The value 7 = 0.4 eV was used 
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in the error metric (Eq. (5)) during the CO, and the spec¬ 
tra were broadened with Gaussians using a full width at 
half maximum (FWHM) of 0.47 eV. The FWHM and 7 
parameters used in the CO were determined by trial and 
error to obtain a suitable balance between the energy and 
intensity sensitivities. The parameters were not specif¬ 
ically optimized and other values that achieve a proper 
balance are expected to yield similar results. For the in¬ 
tegration limits, w m i n = 0 eV was used, and w max was set 
to 5.2 eV, 5.4 eV, and 6.7 eV for Cu 2 , Ag 2 , and Au 2 , re¬ 
spectively. The w max limits were estimated from the grid- 
reference calculations so that spurious box-state transi¬ 
tions are excluded. The same parameters are also used 
in Sec. IV A for calculating the errors with respect to the 
grid reference. The extension phase was terminated when 
the error between consecutive spectra was smaller than 
0.011. The maximum angular momentum in the basis set 
was set to l = 3, allowing for s-, p-, d-, and f-type GTOs 
to be generated by the algorithm. The effect of higher-/ 
basis functions is expected to be insignificant. 

The LCAO spectra of homoatomic dimers were 
calculated within Casida’s linear-response TDDFT 
formalism 63 by including the full eigenstate spectrum 
corresponding to the finite basis set and averaging over 
the longitudinal and transversal components. Other sys¬ 
tems and all grid references were calculated with time- 
propagation TDDFT 62 using a weak 5-pulse and a time 
step of 10 as. The Perdew-Burke-Ernzerhof (PBE) xc 
functional 81,82 in the adiabatic limit was used in all the 
calculations, unless otherwise stated. All the calcula¬ 
tions were done as spin-paired. The following GPAW- 
specific parameters were used. LCAO mode: grid spacing 
h = 0.3 A and minimal vacuum size around the system 
d v ac = 6 A. Grid mode: h = 0.25 A and e/ vac = 8 A. 
The accuracy of the Hartree potential evaluation within 
the simulation cell was improved by employing multipole 
corrections to the potential . 83 The grid mode uses the 
real-space grid for describing the wave functions, which 
explains the smaller h and larger d vac values needed for 
converged results. In the LCAO mode, the numerical ba¬ 
sis functions are described on their specific radial grids 
and the uniform real-space grid is used only for repre¬ 
senting the density and potentials . 66 

In the transferability tests (see Secs. IVB-IVD), the 
spectra were convoluted by a Gaussian broadening with 
FWHM = 0.20 eV. The value 7 = 0.25 eV was used in 
the error metric to obtain a reasonable energy/intensity 
sensitivity. These parameters yield the same error for 
a spectrum shift of 0.1 eV and an intensity change by 
20 % for a test spectrum with a single absorption peak. 
The integration limits in the error metric were set to 
Wmin = 0 eV and w max = 5 eV to probe the visible-near- 
UV region of spectrum. 


IV. RESULTS 
A. Generation of basis sets 

The basis sets were generated independently for Cu, 
Ag, and Au by using the photoabsorption spectrum 
of the corresponding homoatomic metal atom dimer as 
the completeness-optimization target. The dimer bond 
lengths were optimized with GPAW by using the default 
dzp basis sets and the PBE functional. The obtained 
bond lengths for Cu 2 , Ag 2 , and Au 2 are 2.23 A, 2.58 A, 
and 2.55 A, respectively. The used values agree well with 
the experimental values, 2.22 A, 2.53 A, and 2.47 A, 
respectively . 84 

The progression of the CO procedure for the silver 
dimer is illustrated in Fig. 1 for different choices of the 
underlying NAO basis sets. In the extension phase, the 
error in the LCAO calculation decreases rapidly. Once 
the complete basis set (CBS) limit has been reached, the 
least-important primitives are pruned out one by one 
in the reduction phase. During many sequential steps, 
the reduction-phase basis sets yield more accurate re¬ 
sults with less basis functions than the ones from the 
extension phase. The progression of the algorithm is not 
completely monotonic, because the optimized property is 
not variational . 43 

For the rest of the work, we focus on the NAO- 
sz+NGTO basis sets . 85 Then, only a minimal NAO basis 
set is included, so that CO produces all the polarization 
functions necessary for describing the chemical ground- 
state environment as well as the excited state charac- 
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FIG. 1. The progression of the CO procedure for Ag 2 during 
the extension and reduction phases as illustrated by the error 
with respect to the grid-calculated reference spectrum (Eq. 
(5)). Different underlying NAO basis sets are used: 1) “none”: 
only NGTOs, 2) “sz”: single-^ basis of 4d and 5s orbitals, 3) 
“szsvp”: single-^ basis of 4d, 5s, and 5p orbitals, 4) “dzdvp”: 
double-^ basis of 4d, 5s, and 5p orbitals. The cyan circles 
mark the NAO-only basis sets. 
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teristics. Contrary to NAO-only basis sets, the NAO- 
sz+NGTO basis sets are completely general. The occu¬ 
pied orbitals used in the NAO-sz basis can be generated 
for any element, whereas the first unoccupied orbitals 
may not be bound, in which case they cannot be directly 
included. Additionally, the comparison between “sz” and 
“szsvp” reduction series in Fig. 1 indicates that the use 
of the unoccupied p-orbital in the NAO basis does not 
result in a large difference, at least for the silver dimer. 

Furthermore, the NAO-sz+NGTO basis sets are ex¬ 
pected to be numerically efficient. We note from Fig. 1 
that these basis sets indeed perform better than either 
the NAO-dzdvp+NGTO sets or the pure NGTO sets. 
The reasons are that the NAO-dzdvp basis sets contain 
extra NAO functions, which are not free to be optimized 
by the CO method, and that the pure NGTO basis sets 
do not have the advantage of the minimal NAO-sz basis 
sets. However, the NAO basis sets depend on the used xc 
functional whereas the pure NGTO basis sets could be 
more transferable across different functionals. Still, the 
underlying NAOs can be easily changed with the func¬ 
tional. This approach is presented in Sec. IV D. 

To illustrate the generated basis sets, we present two 
NAO-sz+NGTO basis sets for Cu, Ag, and Au in Fig. 2. 
The basis sets with N = 17 are expected to yield decently 
accurate results (see Fig. 1) and the N = 36 / N = 37 
ones are close to the CBS limit. Henceforth, we refer 
to the N = 17 basis sets as “CO-1” and the N = 36 / 
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FIG. 2. NAO-sz+NGTO basis sets generated by the CO al¬ 
gorithm. Basis sets with 17 and 36 or 37 functions for Cu, 
Ag, and Au are shown. 


N = 37 basis sets as “CO-2”. The similarity of the basis 
sets across the studied elements is evident. In the CO-1 
basis sets, the NAO-sz basis is augmented by two diffuse 
p-type NGTOs and a single d-type NGTO. In the CO-2 
case, additional NGTOs are included on all shells. The 
Gaussian exponents of the NGTOs are similar across all 
studied elements, which is expected due to their closely 
related chemical characteristics. Yet, even though the 
shown basis sets are similar for all elements, the whole 
basis set series are not the same, as functions are pruned 
out in different orders in the reduction phase. For ex¬ 
ample, there is no N = 36 basis set for silver. Note also 
that the most diffuse p-type NGTOs are at the constraint 
limit imposed for NGTOs (see Sec. IIIB). 

B. Transferability of basis sets to larger clusters 

The usefulness of basis sets depends on their trans¬ 
ferability to different chemical environments. Here, we 
consider the transferability of the basis sets obtained 
in the previous section to larger systems by using ho- 
moatomic icosahedral clusters of 13, 55, and 147 atoms 
as test cases. The clusters were constructed by adding 
icosahedral Mackay layers one by one around a central 
atom. The structures were relaxed with GPAW by us¬ 
ing the default dzp basis sets and the PBE functional, 
but their icosahedral symmetries were not significantly 
disturbed. 85 

We show in Fig. 3 the error in the photoabsorption 
spectra (Eq. (5)) for the clusters. We observe a nearly 
monotonic increase in the accuracy with increasing basis 
set size. The magnitude of the error is similar between 
different elements, and the error tends to decrease when 
the system size increases. 

In Fig. 3, we also show for comparison the errors of the 
GPAW default dzp basis set and the NAO-dzdvp basis set 
that has been used in a previous study. 25 The default dzp 
basis set is unsuitable for describing the response, which 
is due to its lack of diffuse p-functions. 25 The NAO-dzdvp 
basis set provides an equivalent or better accuracy than 
the NAO-sz+NGTO basis set of similar size. However, in 
contrast to the NAO-only basis sets, the basis sets gen¬ 
erated in the present work allow for further, systematic 
improvements in accuracy beyond that of the NAO-only 
basis sets. 

The insets in Fig. 3 illustrate how the spectra look for 
the 55-atom clusters calculated with basis sets of different 
size. We observe that all the spectral features are mostly 
correct with the basis sets of 17 or more functions per 
atom. The CO-1 spectra suffer from a blue-shift of 0.1 
0.2 eV and the largest improvement in the spectrum when 
growing the size of the basis set comes from a red-shift 
towards the converged spectrum. The CO-2 basis sets of 
Cu and Ag are already at the CBS limit, as the positions 
of the spectral peaks coincide with the grid references 
within 0.05 eV. For Au, the convergence is slow after 
N ss 24, but also there the CO-2 spectrum is near the 
grid reference except for a few details. 
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C. Transferability of basis sets to nanoalloys 

Now, we consider the transferability of the basis sets 
that were optimized for homoatomic dimers to het¬ 
erogeneous metal clusters, where the chemical environ¬ 
ments are different from the homogeneous systems. As 
test systems, we take the core-shell clusters Ag 13 Cu 4 2 , 
Cu 13 Ag 42 , Ag 13 Au 42 , and Au 13 Ag 42 , which consist of 
icosahedral 13-atom cores and single 42-atom Mackay 
layers around the cores, as well as two icosahedral alloys, 
Cu 14 Ag 20 Au 21 and Cu 18 Ag 17 Au 20 , which were generated 
by random substitution of atoms in the 55-atom icosahe¬ 
dral geometry. The clusters were relaxed analogously to 
the homoatomic clusters in Sec. IVB. 85 

The photoabsorption spectra of the alloy clusters are 
shown in Fig. 4. The CO-2 spectra are again in excellent 
agreement with the grid reference. The smaller CO-1 
basis set results in a 0.1-0.2 eV blue-shift of the spectra. 
Due to the lower symmetry and breaking of degeneracies 
the disordered clusters have fewer sharp spectral features 
than the core-shell clusters. This is also reflected in that 
only small differences between the CO-1 and CO-2 spec¬ 
tra can be seen. 

The spectra calculated with the default dzp basis sets 
are also shown in Fig. 4. As in the case of homoatomic 
clusters, these basis sets are unable to sufficiently de¬ 
scribe the photoabsorption spectra. The trends in the 
results obtained with the NAO-dzdvp basis sets (not 
shown) are similar to those in the homoatomic case. 



FIG. 3. Transferability of the generated NAO-sz+NGTO ba¬ 
sis sets to icosahedral homoatomic coinage metal clusters of 
13, 55, and 147 atoms. The off-line filled markers indicate 
results calculated with the GPAW default dzp (N = 15) and 
NAO-dzdvp basis sets (N = 18). The straight gray line is 
drawn to ease the visual comparison. The insets present pho¬ 
toabsorption spectra of 55-atom clusters calculated with basis 
sets of increasing size. The grid reference (thin black line) and 
the LCAO spectrum with the default dzp basis (dotted line) 
are also shown in the insets. 
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FIG. 4. Transferability of the basis sets to alloy clusters. The 
gray shading is applied to the grid-reference spectra. 






























































D. Transferability of basis sets to different xc functionals 


Different xc functionals yield different shapes of the 
KS-orbitals. We studied how the generated NGTOs 
transfer across various xc functionals. The basis sets 
were constructed by augmenting the NAO-sz basis set 
corresponding to the chosen functional with the PBE- 
optimized NGTOs without any further modification. We 
used the Ag 55 cluster as the test system. The results for 
the local density approximation (LDA), 86-88 the Becke- 
Lee-Yang-Parr (BLYP) functional, 89-91 and the solid- 
state modification of the Gritsenko-van Leeuwen-van 
Lenthe-Baerends potential (GLLB-SC) 92,93 are shown in 
Fig. 5, where also the PBE results of Fig. 3 are repeated 
for reference. While the LDA, PBE, and BLYP function¬ 
als predict similar spectra, the GLLB-SC spectrum has a 
distinct shape and stronger intensity. This is due to the 
d-electron screening in coinage metals that is captured 
correctly by the GLLB-SC functional. 25,94,95 The default 
dzp basis sets, shown for comparison in Fig. 5, reproduce 
the strong intensity difference between the GLLB-SC and 
the other functionals but are inadequate to describe the 
detailed structure of the spectra, failing to agree with the 
grid references. In contrast, the CO-1 basis sets mostly 
capture the detailed differences between the xc function¬ 
als, and the accuracy of the CO-1 basis sets is similar 
with all the studied functionals. The CO-2 results are 
again in excellent agreement with the grid references in 
all cases. Altogether, the results illustrate notable trans¬ 
ferability of the PBE-optimized CO basis sets to diverse 
xc functionals. 


E. Computational performance 

A major advantage of LCAO calculations is that their 
computational cost is smaller than that of, e.g., real- 
space grid calculations. However, the advantage de¬ 
creases if large basis sets must be used. To understand 
this important issue, we discuss the effect of basis set 
size on the computational cost. Fig. 6 presents the time- 
propagation run-times of the generated basis sets as cal¬ 
culated for the Ag 55 and Ag 147 clusters with 48 and 96 
processors (cores). 96 Depending on the case, enlarging 
the basis set from CO-1 to CO-2 (i.e., from N = 17 to 
N = 37) increases the computational cost by a factor of 
2 to 5. 

The CO-2 calculations are still 5 to 8 times faster than 
the grid-reference calculations, while the differences be¬ 
tween the results are minimal as observed in the previ¬ 
ous sections. The calculations with the decently-accurate 
CO-1 basis set are 10 to 40 times faster than the cor¬ 
responding grid calculations. 9 ' In the Ag 147 cluster the 
speed-ups are in most cases larger than in the smaller 
Ag 55 system. This is because the studied systems are 
relatively small for the LCAO mode but relatively large 
for the grid mode, i.e., when doubling the number of pro¬ 
cessors from 48 to 96, the grid mode has excellent scal¬ 
ing, whereas for the LCAO mode the benefit from the 
larger number of processors is minor, especially for the 
small Ag 55 . Thus, the speed-up factors are expected to 
be even higher when the system size is further increased 
and the LCAO mode is able to fully take advantage of 
all the available processors. 25 

In addition to the number of basis functions, the com¬ 
putational performance is greatly affected also by the 
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FIG. 5. Transferability of the PBE-optimized basis sets to 
different xc functionals. Results for Ag 55 are shown. The 
GLLB-SC spectrum has been multiplied by a factor of 0.5. 
The gray shading is applied to the grid-reference spectra. 


FIG. 6. Time-propagation run-time of the generated basis set 
series, shown for the Ag 55 and Ag 147 clusters calculated with 
48 and 96 processors (cores). The run-time of the default 
dzp basis set is shown for comparison. The run-times have 
been normalized to the run-times of the corresponding grid 
references calculated with the indicated number of processors. 
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spatial extent of the basis functions. The default dzp 
basis set does not include diffuse functions that are im¬ 
portant for describing the plasmonic response, leading 
to smaller run-time than the CO basis sets of similar 
size. Within the generated basis set series, the effect of 
basis function extent is seen as a staircase-like behav¬ 
ior in Fig. 6. For example, the addition of the f-type 
NGTO (i.e., 7 additional basis functions per atom) in¬ 
creases the computational time only slightly, because the 
added functions are short-ranged. On the other hand, 
the addition of the diffuse s- or p-type NGTO (i.e., 1 or 
3 basis functions per atom, respectively) affects the com¬ 
putational time more clearly due to the functions’ over¬ 
lap with functions on nearby atoms. This effect is pro¬ 
nounced in the larger 147-atom cluster. More aggressive 
truncation of the NGTOs might yield better computa¬ 
tional performance, but to ensure minimal deterioration 
in the accuracy, it may require re-optimization of the ba¬ 
sis sets with the truncation explicitly taken into account 
in the CO routine. 


V. CONCLUSIONS 

In this work, we have addressed the issue of basis set 
completeness in time-dependent density functional the¬ 
ory calculations. We have presented the extension of 
the completeness-optimization paradigm to the genera¬ 
tion of hybrid NAO+NGTO basis sets, and used it to 
parametrize high-accuracy basis sets for nanoplasmonics 
calculations. We have demonstrated the performance of 
the scheme for the coinage metals Cu, Ag, and Au, which 
are experimentally interesting but computationally chal¬ 
lenging due to their semi-core d-electrons that need to be 
modeled in simulations. We have shown that the gener¬ 
ated basis sets are transferable to simulations of various 
metal nanoparticles and nanoalloys as well as to diverse 
xc functionals. 

The results presented in this work are already promis¬ 
ing, but further improvements of the scheme are still pos¬ 
sible. For instance, the error metric used in the present 
work may not be optimal. The metric does not discrim¬ 
inate between different excitations, looking only at the 
aggregate intensity. This may result in spuriously small 
error values due to interference of different excitations. 
The use of an error metric that examines the convergence 
of the excitations one by one might yield even better basis 
set series. 

Another approach deserving further development 
would be to revise the reference systems against which 
the basis sets are optimized. In the present work, ac¬ 
curate and systematically improving basis sets were ob¬ 
tained by optimizing the basis sets for homoatomic metal 
atom dimers. The generated sets were demonstrated to 
be transferable to larger as well as heterogeneous sys¬ 
tems. However, it might be interesting to optimize ba¬ 
sis sets for extended systems, instead. In a dimer, both 
atoms are “on the surface”, which results in the generation 


of diffuse functions to model the exponentially decaying 
density tails. In the solid state there is no exponential 
decay, and diffuse functions are often unnecessary. Never¬ 
theless, our results indicate that the dynamical response 
is sufficiently captured already by the dimer for plasmon- 
ics calculations in larger nanoparticles. 

The main advantage of the LCAO approach is that 
satisfactory results can be obtained much faster than 
with, e.g., grid-type approaches. The main problem of 
the LCAO approach with NAOs (compared to GTO ba¬ 
sis sets of quantum chemistry) has been the scarcity of 
systematically better basis sets. This issue has been ad¬ 
dressed in the present work. 

Although used here for nanoplasmonics, the 
completeness-optimization approach is completely 
general, being applicable to any property at any level of 
theory, also beyond DFT (see, e.g., Refs. 43 and 58). For 
this reason we expect this approach to be widely useful 
in materials modeling by electronic structure methods, 
allowing for large-scale simulations with better control 
on their accuracy. 
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